The classical and quantum dynamics of the inhomogeneous Dicke model and its 

Ehrenfest time 
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We show that in the few-excitation regime the classical and quantum time-evolution of the inho- 
mogeneous Dicke model for N two- level systems coupled to a single boson mode agree for N S> 1. In 
the presence of a single excitation only, the leading term in an 1/iV-expansion of the classical equa- 
tions of motion reproduces the result of the Schrodinger equation. For a small number of excitations, 
the numerical solutions of the classical and quantum problems become equal for N sufficiently large. 
By solving the Schrodinger equation exactly for two excitations and a particular inhomogeneity we 
obtain 1/iV- correct ions which lead to a significant difference between the classical and quantum 
solutions at a new time scale which we identify as an Ehrenferst time, given by tb = y/N/ (g 2 ), 
where \J (g 2 ) is an effective coupling strength between the two-level systems and the boson. 



I. INTRODUCTION 

The recent experimental advances on cold atoms in 
optical cavities 1 , Bose- Einstein condensation of exciton 
polaritons 2 , and observation of vacuum Rabi oscillations 3 
in semiconductor microcavities renewed interest in light- 
matter interaction in the quantum coherent regime. 
These studies were motivated by an observation made 
by Dicke^ long ago who realized that radiation from N 
identical two-level systems (spins 1/2) cannot be treated 
as a sum of N independent radiative processes but rather 
as a collective quantum phenomenon that involves all N 
spins and a photon mode even on the level of perturbation 
theory. Also, several schemes based on light-matter inter- 
action to couple spatially separated spins that had been 
originally proposed as an element of a quantum com- 
puting device^— were recently improved by a suggestion 
to use qubits constructed out of many spins to enhance 
coupling with the optical mode^ due to the superradiant 
effect. 

For instance, considerable attention was paid exper- 
imantally to the -\//V-enhancement of the light- matter 
couplingiiiS. In typical set-ups the spins are spatially 
separated, therefore the excitation energies of different 
spins may be different as they are affected by local forces 
that typically vary across the sample. The coupling 
strength to the light mode also varies as different spins 
are located at different positions of the mode due to a 
different amplitude of the electromagnetic field. Under- 
standing of such inhomogeneities is important to find the 
practical limitations on the decoherence time of the sys- 
tem when, for instance, one designs a quantum comput- 
ing device^'*'— Also, the inhomogeneities are unavoid- 
able and should be important in a system like a semicon- 
ductor quantum dot optical amplifier or laser— 

On the theoretical side, the homogeneous Dicke model, 
which describes a bath of N equivalent spins- 1/2 with 
excitation (Zeeman) energy e coupled to a quantized 
bosonic mode u) with the same coupling constants g, was 
diagonalized exactly in Ref. 15. The influence of inho- 
mogeneities of the coupling constants gj and Zeeman en- 



ergies 6j on the single excitation dynamics was analyzed 
exactly in Refs. EMl It was shown that the boson oc- 
cupation oscillates in time with a single Rabi frequency 
^ = \/N (g 2 ), where yj (g 2 ) is an effective coupling when 
only the coupling constants gj are inhomogeneous but 
with constant Zeeman energies. If the Zeeman ener- 
gies ej are also inhomogeneous but spread narrower than 
the threshold given by f2 this single frequency acquires a 
small Lamb-like shift, whereas for a spread exceeding CI 
the boson decays completely in time. 

In this paper we show that the solution to the classi- 
cal Hamilton equations of motion matches the solution of 
the time-dependent Schrodinger equation when the num- 
ber of spins is large, i.e. N ^> 1, while the number of 
excitations p is still small, i.e. p -C N. For a single ex- 
citation (p — 1) the leading order in an l/A^-expansion 
of the classical equations agrees with the quantum one. 
For a few excitations such correspondence does not hold, 
but for p — 2,3 the numerical solutions of both classi- 
cal equations of motion and Schrodinger equation agree 
for N 3> 1. It is plausible to assume that in leading 
l/A?"-order the same correspondence holds for p > 3. The 
numerical treatment of the Schrodinger equation with a 
large number of spins is possible since the Fock space 
scales only as a power of N (N 2 , iV 3 , . . . ) in the few- 
excitation subspaces. 

As the classical equations of motion for p > 1 can also 
be mapped on the Schrodinger equation in the single ex- 
citation subspace in leading 1 /iV-order the already avail- 
able quantum result can be used to analyze the classical 
equations of motion for few excitations (p -C N) . For p 
excitations with p > 1 we obtain the dynamics by simply 
rescaling the solution derived ir>ii by p. This extends 
the single-excitation quantum solution to the case of few 
excitations when N ^> 1. 

To assess the validity of the classical approximation 
for p > 1 excitations we solve the Schrodinger equa- 
tion exactly in the two-excitation subspace with inho- 
mogeneity in the coupling constants only and compare it 
with the classical solution. When N is small both solu- 
tions are completely different. For large N we perform 
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an 1/A-expansion of the quantum solution and recover 
the classical result in leading order. Subleading 1/N- 
corrections cause deviations between quantum and clas- 
sical dynamics which become significant at a large time 
scale te = yN) (g 2 ) for p <C N. We refer to this time 
scale as an Ehrenfest time, defined here as the time where 
the quantum dynamics starts to differ from the classical 
dynamics. 

Also, having found a quantum solution for p = 2 we 
study it separately and in particular compare it with the 
p = 1 quantum dynamics. We find that inhomogeneity 
of the coupling constants results in a different spectrum 
when N is finite: in the subspace with p = 1 there is 
only one harmonic mode with a single frequency in the 
time-dependent occupation number of the boson, and for 
p = 2 there are N discrete harmonic modes that form 
a continuum spectrum in the limit of large N. Such a 
mechanism can lead to destructive interference, thus to 
decay, of the excitations caused solely by the inhomo- 
geneity of the coupling constants when p > 1. But, as 
pointed out already, for p = 2 we find that the leading 
l/N-tevia recovers the single frequency dynamics in ac- 
cordance with the classical solution. The decay due to 
inhomogeneous coupling constants thus manifests itself 
only in the first subleading 1/iV-correction. We find that 
this contribution is an oscillatory mode with frequency 
|fi and a slowly decaying envelope. The decay behavior 
is essentially non-exponential with a long power-law tail 
and the decay time is r g ~ y/N/ (g 2 ), where y 7 (g 2 ) is a 
characteristic coupling. This decay occurs on the same 
time scale as the Ehrenfest time te defined above. Thus, 
it can be described correctly only by the Schrodinger 
equation (and not by the classical one). 

In our theoretical analysis we assume the following 
ideal experiment. The spin bath is prepared in the 
ground state, e.g. dynamically or by the thermal cool- 
ing. The non-equilibrium dynamics of the boson is then 
initialized by a short radiation pulse from an external 
source which populates the boson mode with a few exci- 
tations like ir*i&~— . The dissipation of the boson mode, 
e.g. leakage of the photons through the mirrors that de- 
fine an optical cavity can be used to detect the dynamics, 
similarly to the measurements performed on semiconduc- 
tor quantum well microcavitie o 2 i 3 i 21 , for the limiting case 
where the cavity leakage time exceeds the internal time 
scale. 

The rest of the paper is organized as follows. In Section 
II we discuss general properties of the inhomogeneous 
Dicke model. In Section III we quote the already known 
solution to the Schrodinger equation in the single exci- 
tation subspace. In Section IV we construct the classi- 
cal analog of the inhomogeneous Dicke model. Section V 
contains the exact solution of the Schrodinger equation in 
the two-excitation subspace for the inhomogeneous cou- 
plings only. In Section VI we compare the numerical so- 
lution of the classical and the quantum equations of mo- 
tions for two and three excitations in the limit of many 
spins. Section VII contains a discussion of applicability 



of the classical approximation. In the Appendix we give 
some details on the calculation of the 1/iV-correction. 

II. INHOMOGENEOUS DICKE MODEL 

The Hamiltonian for the Dicke model that describes 
the interaction between a set of N spins 1/2 with excita- 
tion energies ej and a single bosonic mode of frequency 
to is given by 

JV N 

H = uitfb + ejS! + £ 9i {S+b + Srtf) , (1) 

3=1 3=1 

where Sf = SJ ± iS? , Sj are spin 1/2 operators, b mM 
are the standard Bose annihilation (creation) operators. 
The coupling constants gj are typically given as dipole 
matrix elements and thus are, in general, complex num- 
bers. Since their phases can be eliminated by a unitary 
transformation, we treat gj as real and positive numbers. 

In the present paper we assume that the boson mode 
is tuned in resonance with the spins (ej) — w, where 
(...) = Y^j ■ ■ ■ /N ■ If the boson mode is strongly de- 
tuned, \(ej) — (*}\ ^> y (.g 2 ), the interaction between 
them is weak and the model Eq. ([1} can be analyzed 
perturbativelyii. Also note that the inhomogeneities of 
gj and/or ej forbids to represent the Hamiltonian Eq. 
(Q]) in terms of the total angular momentum operators 
J a = Ej S f, at = x,y,z. . 

The total number of spin-boson excitations, L = n + 
Ylj Sj, is conserved by the model Eq. (JXJ) , where n = b'b 
is the bosonic occupation number. The eigenvalue c of L 
labels the subspace of the Hamiltonian with a given total 
number of excitations. 

We restrict ourselves to a small number of excitations, 
p -C N. In the following we assume that the spins can 
be prepared in the ground state with each spin in its low 
Zeeman state. The bosonic mode is assumed to be occu- 
pied by p bosons initially, the time evolution is restricted 
to the subspace with c = —N/2 + p. Then the leakage 
of the boson mode to the outside world can be used to 
monitor the time dynamics of the system by detecting 
the leaked mode at given subsequent instances in time. 

III. SINGLE EXCITATION 

The time dynamics of Eq. (CD) for a single excitation 
was analyzed in detail in Ref. 17. Here we only quote the 
explicit form of the corresponding Schrodinger equation 
and the main results derived from it. 

The time evolution is restricted the the subspace with 
c = —N/2 + 1 and is described by the general state 

N 

= a (*) |4, + |^,0>, (2) 

3=1 
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where a (t) and /3j (t) are normalized amplitudes, 

\ a W| 2 + J2j \Pj {t)\ 2 = 1) 01 finding either a state with 
one boson and no spin excitations present or a state with 
no boson and the j* -spin excited (flipped). As initial 
condition we will assume throughout (with one excep- 
tion discussed at the end) that initially only bosonic ex- 
citations are present while each spin is in its individual 
ground state, i.e. a (t = 0) = 1. The state from 
Eq. |2} describes then the time evolution of an initial 
product state 1) into an entangled state formed by a 
coherent superposition of N+l states, where each |JJ-tj) 0) 
contains an excited spin and no boson. This entangled 
state can be viewed as a (para-) magnon state in the 
uniform limit. In other words, the initial bosonic excita- 
tion gets coherently spread out over the entire system in 
course of time. 

Inserting \*f? (t)} from Eq. into the time-dependent 
Schrodinger equation we get 

-ia{t) = X>&(*)> ( 3 ) 

3 

-0 k (t) = (ej -Lj)p k (t)+g k a(t). 

This set of coupled equations can be solved explicitly 
via Laplace transformation. We use the same approach 
to solve the Schrodinger equation in the two excitation 
subspace in Section V of this paper. 

If the number of spins is large, N ^> 1, the sum over 
j in the exact solution of Eq. ([3]) can be substituted by 
an integral. In this continuum limit the discrete set of 
6j and gj become continuous variables characterized by 
distribution functions Q (g) and -P(e). Any distribution 
function of g results only in a renormalized coupling con- 
stant \J {g 2 } and the dynamics of the boson is not affected 
in any other way. 

Different distribution functions of e result in qualita- 
tively different regimes of the dynamics. Let us choose 
P (e) as a rectangular pulse shape of width A centered 
around w, 

P(e) =6>(-e + cj + A/2)(9(e-u; + A/2), (4) 

where 9(x) is the Heaviside step-function. It was shown 
that when the inhomogeneit y is be low a certain thresh- 
old, A/fl <C 1, where £1 = J N (g 2 ) is the collective Rabi 
frequency, the boson excitation, (n) = \a(t)\ 2 , does not 
decay, i. e. 

(n(t)) = cos 2 (fit) . (5) 

The corrections to this result are small and are on the 
order of A/SI. In the opposite limit, A/fl ^> 1, the spins 
act as the thermal bath at zero temperature. The bosonic 
excitation decays completely and exponentially, 

(n(t)> = exp (-t/ta) , (6) 

with the decay time t% = 2A/irVl 2 . In the intermediate 
regime, A ~ fi, the decay is partial and the decay law is 
a combination of exponential and inverse-power laws. 



IV. CLASSICAL ANALOGY 

Here we construct a classical version of the inhomoge- 
neous Dicke model. Using Dirac's analogy^ we change 
the boson operator b in the model Eq. (JXJ) to a classical 
complex variable a = a x + ia y , and the spin operators 
Sj to a set of N vectors Cj = (C?,CJ,C|) of length 
|Cj| = 1/2. These classical degrees of freedom obey the 
Poisson bracket relations which are obtained from the 
bosonic and spin commutation relations via the ansatz 
[j] — * Del- [Ca,Cp] cl = -e a/ 3 7 C 7 , and [a,a*] cl = i. 

The Hamilton equation of motion for the j th spin, 
Cj = [H, Cj] ,, is a Bloch equation 

Cj = 1* x Cj, (7) 

where the in-plane component of the effective magnetic 
field is the complex bosonic field, and the perpendicu- 
lar component is the single spin excitation energy, = 
(2gja x , 2gjdy, tj — uj). The Hamilton equation of motion 
for a is a feedback to the bosonic field from the in-plane 
component of all spins, 

3 

where CJ = CJ — iCj . Generally, these differential equa- 
tions can be solved numerically with the initial conditions 
Cj (0) = (0, 0, -1/2) and a (0) = ^/p to obtain the time- 
dependent solution Cj (t) and a (t) explicitly. The time- 
dependent value of the bosonic field is n (t) — \a(t)\ 2 . 

For a small number of excitations, p <C N, Eq. ([7]) 
simplifies. The quantity L = \a\ = ~ N/2 +p is 

conserved during the evolution governed by Eqs. ([71 [5J. 
Thus, at any instance of time J2j Cj ~ —N/2, i. e. if 
the dynamics starts with only a few bosonic excitations, 
the spins cannot 'flip' during the evolution. Using the 
approximation Cj (t) ps —1/2, the equation for CJ (t) 
drops out from Eq. (J7J) and the remaining two equations 
are 

Cj = -i (ej - uj) Cj ~ ig 3 a. (9) 

When p = 1, the Hamilton Eqs. (0 [9|) with the initial 
conditions above coincide formally with the Schrodinger 
Eq. ([5]) . By direct comparison we can establish the corre- 
spondence between the quantum mechanical amplitudes 
and the classical variables: the classical field a is the am- 
plitude a and the in-plane component of the spin vector 
Cj is the amplitude j3j . Note that the classical spins are 
not averages of the spin operators, (Sj) = 0, but instead 
they are connected with the quantum mechanical ampli- 
tudes. The solution of the dynamical Eqs. (JH1 [HJ) is the 
same as the solution of Eq. ([3]) . 

When the number of excitations at the initial time is 
p > 1, but still p <C N, the mapping of the classical 
equations on Eq. ([3]), approximation (JSJ), still holds but 
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the initial condition for Eq. is different: a (0) = yfp. 
This difference results then only in a renormalization of 
(n(t)) [obtained from Eq. by p. 

When the number of excitations is large p > N, the 
z-components of the classical spins deviate significantly 
from their initial values during the evolution and the ap- 
proximation j9j) is not valid. In this regime the classical 
Eqs. (0 [H} have to be solved numerically. 



V. TWO-EXCITATION REGIME AND 
INHOMOGENEOUS COUPLING CONSTANTS 

In this section we consider the dynamics of two exci- 
tations for a system with inhomogeneous coupling con- 
stants gj but constant Zeeman energies ej = co. 

The time evolution of the two excitations is restricted 
to the subspace with c = —N/2 + 2 and is described by 
the general state 

|* (t)> = a (t) |2,4) + ft Mti) + E 7« Mtitj) , 

3 i>3 

( 10 ) 

where a(t), (t), and jij (t) are the normalized ampli- 
tudes, |a(t)| 2 +Ej I ft (*)| 2 +Ei>j |7i/ = of the state 
with two bosonic excitations, a state with one bosonic 
excitation and the j spin excited, and a state with no 
bosonic excitation and the i th and j th spins excited (with 
i =/= j). The amplitude 7^ is defined such that 7^ = if 
j > i- 

The conservation law can be used to simplify the 
Hamiltonian Eq. (p}. We subtract ojL from Eq. JT]), 
which only changes an irrelevant overall phase of |^ (t)}, 
to eliminate the first two terms. Note that the second 
term will not be zero away from the resonance u 7^ e. In- 
serting \*f? (t)) into the time-dependent Schrodinger equa- 
tion we then obtain, 



A. General solution 

We use the Laplace transform, A (s) = J °° dtA (t) e st , 
to solve the set of equations, Eq. (fTTj) . In the Laplace 
domain Eq. (|lip is a set of linear algebraic equations, 

-i(aa-l) = \2V i!(i 3. 

-isfik = V2g k a + Y.j<k 9jlkj + Y,j>k 9j7jk 
-isjkl = (Sfcft + gifik) (1 - hi) , 

that can be explicitly solved. The substitution of a (s) 
and 7fc; (s) as functions of /?/. (s), that are obtained from 
the first and the last lines, into the middle line gives the 
following set of equations for /3k only, 




3X>&-tV2. (14) 



Each /3fe (s) is easily found from the the above equation 
since J^j 9jPj ( s ) on t ne right hand side is the same in 
each equation for all /3& (s). Then the sum is found self- 
consistently and we obtain the solution for the amplitude 



as 



-i\[2g k 



2 5 



l- N (9 2 )l- 3 (- 



gJN 

■ S 2+2g?-N{g*) 



(15) 

where the average is the sum over all spins (...) = 
(y^,j ■ ■ ■ J /N . The other two amplitudes are found from 
the first and the third lines of Eq. (|13j) by substitution 
of the above solution for (3k (s), 



1 / 9\N 

1 \- s 2 +2 g*-N(g*) 

a (s) = 

s 1 -3 



- s 2 +2 g*-N{g2) 



(16) 



-id(t) = \/2^>ft(i), (11) 

-i/3 fc (*) = V2g k a (t) + ^ ^7^ (t) + ^ 5j7jfe (i) , 

j<k j>k 

-*7«(*) = CfffcA (*)+SlM*)) (!-*«)• 

The initial condition, a (0) = 1 and ft (0) = 7^ (0) = 
0, which we further assume corresponds to the doubly 
occupied boson mode at the initial time. The physical 
observable of interest is the time-dependent value of the 
boson occupation number (n (t)}, which can be expressed 
in terms of the amplitudes a (t) and ft (t) as 

(n(t))=2\a(t)\ 2 +J2\fr(t)\ 2 , (12) 

3 

where (...) = (\& (t)|. . it)) is the time-dependent ex- 
pectation value. 



lki(s) = -{g k Pi + giPk) (X - Ski) ■ (17) 
s 

The main focus of our interest will be on Eqs. (If 51 
[TH]) as the observable quantity (n (i)) depends only on 
a (i) and /3& (<) . These time-dependent amplitudes can 
be obtained from Eqs. (|15[ ITS)) by the inverse Laplace 
transform. The analytic structure of Eqs. (fT51 HT))) is 
governed, in general, by a set of poles given by the roots 
of denominators which depend on a particular set of gj . 
For instance, if the number of spins N is small there 
are 2N conjugated complex roots. The inverse Laplace 
transforms of a (s) and /3k (s) will be a sum of N dis- 
crete harmonic modes in contrast to the single excitation 
dynamics where in such a setup there is just a single 
pair of roots independent of the particular set of gj , see 
Section III, and there is only a single harmonic mode in 
the dynamics of the boson occupation number, see Eq. 
([3]). Such a result marks a qualitative difference in the 
dynamics of the single- and two-excitation subspaces. 
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B. Time-evolution in the continuum limit of many 

spins 



go ~ £ to a maximum value g = go, 



In this section we study the limit of many spins, i.e. 
N ^> 1. The sum over j in Eqs. (fTS! \W\\ can be sub- 
stituted by an integral over a distribution function of g, 
Y]- • - ■ N / °° dgQ (g) ■ ■ ■ In the continuum limit some 
poles can merge together, forming branch cuts, and some 
poles can separate themselves from the others. The in- 
verse Laplace transform of the branch cuts will become a 
decay function in the time domain and the separate poles 
will contribute a set of harmonic modes. 

The analytic structure of a (s) and (s) explicitly de- 
pends on the particular form of Q (g). 



Uniform distribution function 



Qi (9) = e(-g+go)e(g-9o + /e 



(18) 



The coupling constants cannot be negative, so £ can vary 
from £ = (e. g. all couplings are the same and are 
equal to go) to £ = go (e. g. the couplings are evenly 
distributed from to go), see Fig. QJi. A useful property 
of this distribution function Qi is that a small and a 
large inhomogeneity can be analyzed on the same footing. 
Another distribution function will be considered in the 
next subsection. 



To be specific we consider a set of coupling constants 
which are uniformly distributed from a minimum value 



Turning the sum in Eqs. (|15l I16|) into an integral and 
using Qi (<?), we obtain 



J 



g]N 



s 2 + 2g 2 - N (g 2 ) 




N(g 2 



■ arctan 



N(g 2 



+ N (.g 2 ) - 2g (go - 



r 



(19) 



Qi(g) 



a) 





■< >• 











9o 



go 




Figure 1: Distribution functions of g that are used to evaluate 
the sums in Eqs. (|15II16[) . (a) The uniform distribution func- 
tion Qi (g) has a maximum coupling strength go and a width £ 
which can vary from (i.e. homogeneous coupling constants) 
to go (i.e. maximally inhomogeneous coupling constants), (b) 
The sawtooth distribution function Q2 (g) describes a non- 
uniform spread of the coupling constants gj from to go ■ 



where (g 2 ) = g 2 - £g + £ 2 A 

The analytic structure of a (s) from Eq. (fll)]) with the 
sum from Eq. (|19[) is the following. There are three 
poles and two branch cuts, see Fig. [5J Thus, the in- 
verse Laplace transform has two contributions a (t) = 
a p (t) + a c (t). One pole is at s — and two poles are 
at s = Usq, where sq = . These are given 

by zeroes of the denominator of the second term in the 



product in Eq. ([Trj)) . Note that so was obtained using a 
1/iV-expansion and is independent of £ in leading order. 
In the first subleading 1/iV-order so depends on £, 



so 



ldN{g 2 )--gl 



when £ = go, and 



s = 2dN(g 2 ) 



-,gl 



(20) 



(21) 



when £ = 0. The inverse Laplace transform of the func- 
tions with poles is a sum over the corresponding residues, 

a p(t) = J2s=o,±is Res » a ( s ) eS *> and ^ g ives 



a p (t) 



1 



cos sot 



1 

N 



8a 



(22) 



Here, the leading term is independent of £ unlike the 



first l/7V-correction, 5a (go) 
fo(0)- 11 



27 
20 



(1 — cos (s t)) and 



i I cos (sot)) . We refer to Appendix A for 
the calculation. 



The expression in Eq. (|19[) has four branch points. 
Two are given by the square root, s — ±isi where 
• s i = V N (g 2 )- The remaining two are given by arctan. 
Solving the equation 



N(g 2 



+ N (g 2 ) - 2g (go - 



±i 



(23) 
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a(s) 



Im(s) 
is 

isi 
is 2 

4 o 

-is 2 
— is i 

-is 



Ms) 



Im(s) 
; is 

is 2 

-is 2 
— is± 

-is 



Figure 2: Analytic structure of the time-dependent quantum 
mechanical amplitudes for p = 2 excitations in the Laplace 
domain, Eqs. (|15l I16[) , in the continuum approximation cal- 
culated using the distribution functions Q\ (g) and Q2 (<?)■ 
Separated dots, 0, ±iso, are poles and the dots, iisi^, con- 
nected by bold lines, are branch points. The bold lines are 
the corresponding branch cuts. 



we 



them as 



= ±is2 where S2 = 



find 

^N(gi)- e ~ 2 5 o (50 - ~ W^ 2 +4.90 (flo-0- 
The first branch cut is chosen as a straight line between 
isi and IS2 and the second branch cut as a straight line 
between — is% and — isi, see Fig. [2] 

The contribution to the inverse Laplace transform from 
the branch cuts is a function of £. When £ = Eq. (fTU)) 
has no branch points. In the £ — > limit the arctan can 
be expanded in the small parameter, then the leading 
term is non-zero and contains no multivalued functions. 
All the higher order terms are proportional to £ and are 
zero when £ = 0. We obtain in this limit a c (t) = 0. 

When £ = go, the integral enclosing the branch cuts, 



(*) = 



32 
3iV2 



X 2 cos 



dx- 



- 1 



l. 2 X ) 



(24) 



contributes only to the second subleading 1/iV-order of 
a (t). Thus, a c (t) is beyond the accuracy of the present 
calculation for all values of £ and it can be neglected 
compared to the leading order correction in Eq. (f2"2"j) . 



The analytic structure of (3 k (s) in Eq. (fT5)) is the same 
as a (s) except that there is no pole at s — 0, see Fig. [5] 
Thus the inverse Laplace transform also has two contri- 
butions, (3 k (t) = (3 p k \t)+(3 c k (t), when £ > 0. One is given 



by the sum over just two residues, s 



three 



= ±iso, instead of 
and yields 



-ig k sin (apt) 



N 



(25) 



where, similarly to Eq. ((2"2"|) . only the first 1/N correction 
depends on £ but the leading term does not, 5 (3 k (go) = 
-igk (2 (g k /gof - 3) sin (s t) /^2N (g 2 ) and 8(3 k (0) = 
jsin(soi) /1\j2N , see Appendix A for the calculation. 

When £ = 0, the branch cuts disappear, (3% (t) = 0, 
similarly to a c (t). When £, = go the analysis of the 
branch cuts is a bit different from above for a c (s). There 
is a singularity in Eq. (fT5|) at s = ±i N (g 2 ) — 2g k , orig- 
inating from the first term in the product in Eq. (1151) . 
which overlaps with the branch cuts. It present a diffi- 
culty if we apply continuum approximation to the dis- 
crete form of (3k (s) in the same way as we did to a (s). 
Cancellation of this singularity by a zero in the denomina- 
tor of the second term in the product in the original dis- 
crete form, Eq. (|15j) . simplifies the analysis. The analytic 
structure of (3k (s) in the continuum approximation does 
not alter. The only difference is a small 1 /TV-correction 
to Eq. (fl"9f . Repeating the same steps as between Eq. 
(fT9| and Eq. (|24f we obtain the following expression for 
the integral enclosing the branch cuts, 



& (*) = " 



2v / 2^5o 2 



2i sin 



dx- 



N< 



2g^x 2 t 



3N 



(26) 



Here, the integral can be simplified by performing an 2go/3N. As a result, the leading 1/N -term is 
1/iV-expansion. This approximation is valid for the ma- 
jority of gk 3> 2g /3N except for a small set of gk <C /o A' 2 
2g /3N, where the maximum value of \(3 c k (t)\ < ^ P c k (*) = - — — j=±= I (t) , (27) 
is small as 1/s/N compared to the majority of gk ^> 



where the dimensionless integral / (t) describes the time- 
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Figure 3: Decay correction to the dynamics of the boson Eq. 
(|28[) for the parameter yj N (g 2 )/go = 5. The solid line is a 
numerical evaluation of the integral, where t 9 marks the time 
scale of the initial decay. The dashed line is the long time 
asymptote Eq. ([30l . 



decay, 



I(t) 



x 2 sin 



dx- 



2g 2 x 2 t 



(fl n (l-f 



- 1 



+ (!*) 



(28) 



and is independent of k. In contrast to a c (t), /3% (t) does 
contribute to the first subleading 1/N-oidei of (3k it) and 
we will analyze it below. 

The argument of the sine in Eq. (|28jl can be expanded 

The 



1/N, y/N(g 2 )-2g 2 x 2 « y/W^) - 



2 2 



leading term yj N (g 2 )t, which is a fast oscillating func- 
tion, can be taken outside of the integral. The second 
term gives a slow decay envelope. This term leads to 
a significantly decay when g 2 t/ y/ N (g 2 ) = 1. Thus, we 
estimate the decay time as 



= VW), 



(29) 



different by a factor of 3, see Fig. [3J as we neglected the 
logarithmic singularity at x = 1 in the initial integral. 

In the sum of two amplitudes a (t) and /3k it) the decay 
shows up only in the fist subleading 1 /TV-order when the 
number of spins is large. The particular form of the decay 
function is rather involved and is not displayed here. 

The decay on the short time scale t < r g is essen- 
tially non-exponential, see Fig. [3J The particular shape 
depends on the particular set of the strongest coupling 
constants gj. However, an estimate of the time scale 
T g ~ \/Nj (g 2 ) is independent of Q (g), as it is based on 
an 1 /TV-expansion only, i. e. the distance between the 
branch points in Fig. [3] is smaller by 1/N compared to 
the distances between the branch points and the poles. 
The power law tail exists due to a bound on the smallest 
gj. The power and the numerical prefactor in Eq. (|30[) 
depend on a particular Q ig), especially on the distribu- 
tion of the smallest gj's as they are responsible for the 
long time behavior. 

The time-dependent occupation number of the boson 
can also be expanded into a 1/iV-series. The leading term 
depends on £ only through the effective coupling yj (g 2 ). 



(nit)) = 2 cos 2 (f^j+ife it). (31) 

The leading l/7V-correction, 6n(t) — Aa (i) 5a it) + 
2 Pj (t) Sf3j it), is qualitatively different for £ = and 
£ > 0. When the coupling constants are homogeneous, 

£ = o, 

l 



5n it) = - ( 2 (1 + cos (s *)) + g (1 



cos(2s t)) (32) 



is an oscillatory function where the second harmonic with 
the doubled frequency 2sq appears in addition to the 
main frequency so of the leading term. For the case of 
maximally inhomogeneous coupling constants, t; = go- 
the function 



see Fig. H At a large time t > y/N/ (g 2 ), I it) has 
a power-law tail. Due to fast oscillations of the sine 
the main contribution to the integral comes from x <C 
{/ N (g 2 ) I \J g 2 t, thus, the spectral function can be ap- 
proximated as x 2 1 ^| In (izf) - l) + (f %)' 

Then, the integral in Eq. (f2"5|) evaluates in terms of an 
error function which we expand, again, in a Taylor series 
in powers of yj N (g 2 ) / g 2 t <C 1, and we obtain 



/(*) = 



X 
2g a t 



cos 



x - ^ ) V 



1 

+ 2 



TTX_ 

2g t 



(sin ixg Q t) + cos ixgot)) , (30) 



where \ — \J N (g 2 ) /go- The shape of the asymptote 
is in qualitative agreement with the explicit numerical 
evaluation of Eq. (1281) . however, the overall amplitude is 



Sn it) = 1.8 sin 2 (s t) - 2 sin (s t) I it) (33) 

contains a decaying contribution, where / (t) is the decay 
function from Eq. (f2"5|) . Up to the time scale r g the term 
sin (sot) / it) ~ sin [sgt) cos is it) is a harmonic mode with 
the third frequency So + Si in addition to Sq and 2s . This 
mode can only be observed if the short-time regime with 
t < T g is accessible to measurement. 



D. Sawtooth distribution function Q2 ig) 

Here we study another distribution of the coupling con- 
stants. Assuming that there are more spins at the nodes 
of the cavity mode, so that the stronger coupling con- 
stants are more favorable, we consider a sawtooth-like 
distribution function with its maximum at the largest 
coupling strength g , Q 2 ig) = 2g/g 2 6 (-g + g ) 9 ig), see 
Fig. 



Replacing the sums in Eqs. (fT51 [To]) by integrals, ■ • • N J °° dgQ (g) . . . , and using Qi (g) we get 



g]N 



s 2 + 2g 2 - N (g 2 ) 



N 



2^ 



+ N(g 2 ) 



+ N(g 2 )-2g 2 



(34) 



r 



where (g 2 ) = g$/2. 

The analytic structure of Eqs. (fl~5| [T6|) with the sum 
from the above equation is the same as with the sum 
from Eq. (flU)) obtained using Q\ (g). There are three 
poles (a (s) has three poles and (3k (s) has only two as in 
the previous subsection) and two branch cuts, see Fig. [2] 
Thus, the inverse Laplace transforms of a (s) and /3k (s) 
also have two contributions, i.e. a (t) — a p (t) + a c (t) 
and /3k (t) = 0u if) + {31 (t) . Similarly to the previous 
subsection, there is a pole at s — and there are two poles 
at s = ±iso, where sq = 2yjN (g 2 ) agrees in leading 1/iV"- 
order with what was obtained in the previous subsection. 
The four branch points, which are due to the logarithm 
in Eq. |34|) . are found from 



+ N(g 2 ) 



+ N (g 2 ) - 2g 2 



0, oo, 



(35) 



S-2 



as s = Hsi t 2 where s% = yj N (g 2 ) and 
y/N (g 2 ) - 2g 2 . These also agree with what we have al- 
ready found in the previous subsection when the coupling 
constants were maximally inhomogeneous, i.e. £ = go. 

The sums over the residues give the main contribu- 
tion to the inverse Laplace transforms. The leading 1/N- 
terms in a (t) and j3k (t) agree with the leading terms in 
Eq. (1221) and Eq. (p?5|) , where (g 2 ) has to be calculated 
using the sawtooth distribution function Q2 (g) instead 
of Qi (g). The contributions from the branch cuts also 
appear only in the first subleading 1/N order. The main 
features of the time decay are similar to that of Eq. . 
Indeed, the decay time T g and the frequency of the fast 
oscillating term in Eq. (|28[) for times t < T g result from 
the same branch points, ±isx,2, as in the previous sub- 
section. 



As the amplitudes a (t) and j3k (t) are similar to the 
ones in the previous subsection, the boson number (n(t)} 
is also given by Eq. (|3"Tj) . The leading l/iV-term depends 
on Q2 (g) only through the effective coupling constant 
yj (g 2 ), and the leading 1 /iV-correction contains a decay 
term. 



VI. SUBSPACE OF TWO AND THREE 
EXCITATIONS 

In this section we compare numerically the solution of 
the Schrodinger equation with the one of the classical 
equations of motion Eqs. ([71 H]) for p — 2, 3 excitations. 
We start from writing down the Schrodinger equation in 
these two subspaces explicitly. 

The time evolution of two excitations is restricted to 
the subspace with c = —N/2 + 2 and is described by the 
general state Eq. (jTUJ) . The Schrodinger equation for an 
arbitrary set of ej and gj in this subspace is similar to 
Eq. (HU, 



(36) 



-0k = (ek-u)Pk + V2gka + '^2gjlkj+'^2gj'Yjk 

j<k j>k 

-ijki = [(e* - w) + (ej - u)]"fki + {gkPi + giPk) (1 - Ski) 

with the initial condition a (0) = 1, /3& (0) = 0, and 
lu (0) = 0. 

The subspace of three excitations is labeled by c = 
—N/2 + 3 and is described by the general state, 



J 



|*(t)) =a(t)|3^>+X>|2^ti}+X>«l 1 >^it J -)+ E ^iltUtrtitj), 



i>j>r 



(37) 



;lh 



where a (t) , fij (t) , 7^ (t) , and rj r ij (t) are the normal- spin excited, a state with one bosonic excitation and the 
ized amplitudes, |a(i)| 2 + £\ (t)\ 2 + Y Jl>] + 
^2i>j>7- IVrij] — 1) °f the state with three bosonic exci- 
tations, a state with two bosonic excitations and the j th 



and j spins excited, and a state with no bosonic ex- 
citations and the i th , j th and r th spins excited. The am- 
plitude jij is defined such that 7^ = if j > i, and rj ri j 
is defined such that 7? r y = if the inequality i > j > r is 
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not satisfied. The Schrodinger equation in this subspace 
I 



— ia = 

-0k = 

-ijki = 

^Vklm — 



with the initial conditions a (0) = 1, f3 k (0) = 0, and 
lu (0) = 0. The physical observable of interest is again 
the time-dependent boson number which can be ex- 
pressed in terms of the amplitudes a (t), j3j (t), and 7^ (t) 
as 

(n(t)) = 3 \a (t)\ 2 + 2 £ |# (t)\ 2 + £ (t)\ 2 . (39) 
3 3 

Unlike before for p = 1, the classical equations of mo- 
tion Eqs. (J7J H} and the Schrodinger Eqs. (US [35]) are 
not equivalent in the subspaces of p = 2, 3. Solving these 
equations numerically, we compare the time-dependent 
boson number, given by Eqs. (fT2l |3"9")) . for p = 2, 3 with 
the square modulus of the classical field a (t) obtained 
from Eqs. (0 [5J. In the large- N limit the solutions of 
both classical and quantum equations have exactly the 
same form. 

In Fig. Q] we give a detailed comparison of the differ- 
ent types of inhomogeneities characterized by the distri- 
bution functions P (e) and Q\ (g) for a large number of 
spins, i.e. TV = 200. The solutions of Eq. {35) and Eqs. 
(UJ [SJ with p = 2 are compared in Fig. 0] a)-c) , and the 
solutions of Eq. (|38[) and Eqs. ([7J [S]) with p = 3 are com- 
pared in Fig. [5]d)-f). The regime with inhomogeneous 
coupling constants only, characterized by Qi (g) from Eq. 
(fT8f with £ = 50, and P(e) from Eq. g]) with A = 0, 
is shown in Fig. H] b) and e). The regime of strongly 
inhomogeneous Zeeman energies, A/fi = 10, is shown in 
Fig. 0] c) and f ) . Finally, an intermediate regime with 
A/fi = 2.2 is shown in Fig. [4] a) and d). 

As the classical and quantum solutions coincide in all 
regimes, the classical equations can be used to find the 
time-dynamics of the boson occupation number. This is 
quite remarkable since the classical equations are signif- 
icantly simpler to solve than the Schrodinger equation, 
both analytically and numerically. When the number of 
excitation is small, p -C N, the classical equations can be 
mapped to the Schrodinger equation in the one-excitation 
subspace p — 1 in leading 1/iV-order, see Eqs. (JHHU)- The 
Schrodinger equation for this case was already solved. A 
larger number p of excitations changes only the initial 



is 



condition of Eq. © from a (0) = 1 to a (0) = y/p. This 
difference can be accounted for by a simple rescaling by 
p of the bosonic occupation number (n (t)} that was ob- 
tained in the single-excitation subspace for all regimes. 
The explicit results were given in Sec. Ill, see Eqs. (0 
O, but we do not show them in Fig. 01 

In conclusion, the analysis in Ref. [I?] is also applicable 
to the case with more than one excitation, provided p <C 
N and N > 1. 



VII. APPLICABILITY OF THE CLASSICAL 
APPROXIMATION 

The classical approximation in the few-excitation sec- 
tor is exact when the number of spins N is infinite. For 
finite but still large N's, i.e. N ^ 1, the time-evolution 
of the classical system deviates from the quantum system 
by a small amount. The goal of this section is to analyze 
these finite-size deviations 2 ^. 



A. Ehrenfest time 

One way to quantify the difference between the classi- 
cal and the quantum solution is to identify the Ehren- 
fest time te at which they deviate significantly from 
each other. To this end, we compare the boson num- 
ber (n(t)) , obtained from the classical equations (in the 
few-excitation approximation), Eqs. (JHHU), with the ex- 
act quantum solution in the two-excitation subspace ob- 
tained in Sec. V. The solution to the classical equa- 
tions Eq. ([5]) in this regime is n — 2 cos 2 (fit), where 
Q = ^JN (g 2 ) . From Eq. j3l|) the solution to the 
Schrodinger equation is (n) — 2 cos 2 (s^t / '2) in leading 
1/N-oidei. Both solutions are single harmonic modes 
with frequencies that also match in leading 1/A^-order, 
where sq = 2-JN (g 2 ). In first subleading 1/W-order, the 
correction to sq depends explicitly on Qi (<?), see Eqs. 
([HI EH)- For £ = the expansion of Eq. (JUJ) gives 
s = 2^/N{g 2 ) - yfl (g 2 ) /N, and for £ = g expanding 



yfi^gjft, (38) 

3 

(e k - w) /3 k + V?>g k a + V2^ 9jlkj + ^^djljk, 

j<k j>k 

[(e k - u) + (e t - w)]7jw + V2(g k Pi + gi/3 k ) (1 - S k i) + 9jVjki + ^2 9jVkji + ^2 9jVkij, 

j>k>l k>j>l k>l>j 

[(e k + (ei + (e m - u)] Jkl + gkjlm + gilkm + 9mllm 

I 




Figure 4: Numerical solutions of the Schrodinger Eqs. (|36l 138ft and classical equations of motion Eqs. ([7l [Sj> Ior a large number 
of spins, N = 200, and p = 2, 3 number of excitations. The inhomogeneities are characterized by the distributions P (e) and 
Qi (g). In plots a) and d): A/fi = 2.2 and £ = 0; in plots b) and e): A = and £ = go; in plots c) and f): A/fi = 10 and £ = 0. 
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Eq. (gnjl we get s = 2^N (g 2 ) - 3^3 (g 2 ) /ION. Thus, 
at the time scale 

/ N , x 

the phase difference between the two harmonic modes, 
with frequencies Q and sq/2, is comparable to 2tt. Hence, 
the difference between classical and quantum solutions is 
significant for any value of £ at this time scale te, to 
which we refer as Ehrenfest time. The numerical pref- 
actor depends on the particular type of inhomogeneity. 
The explicit calculations in Sec. V for typical distribu- 
tion functions show that these numerical prefactors are 
of order one. 

There is also a 1/iV-correction to the amplitude of 
(n(t)) coming from the Schrodinger equation, see Eq. 
(|33p . This correction contains a decaying contribution 
(proportional to I(t) given in Eq. ([2"5)0 which, again, 
marks a qualitative difference between the quantum and 
classical time-dynamics in the two-excitation subspace. 
In particular, I(t) decays at the characteristic time scale 
T g given in Eq. (|29j) which is equal to the Ehrenfest time 
te introduced above. Thus, we see that this difference 
in the amplitude (although it is only a 1/N-correction) 
is another manifestation of the quantum nature of the 
system where the time-dynamics for times exceeding the 
Ehrenfest time can be described correctly only by the 
Schrodinger equation (and not by the classical one). 

So far we have been using the approximate classical 
Eqs. ([5J O for few excitations, thereby neglecting the 
deviations of the z-component of the classical spins from 
Cj = —1/2. To estimate the quality of this approxi- 
mation, we use the result obtained in Ref. |26| for the 
homogeneous classical system. The solution of the un- 
approximated Eqs. (0 [5]) with gj = go and ej = lj is 
an elliptic function of time2£. When the number of ex- 
citations is small, p -C N, this elliptic function can be 
expanded into a harmonic series with a leading term that 
reproduces the solution of the approximate Eqs. (0 [H]). 
The frequency of the leading harmonic term matches the 
frequency Q — \/N go from Eq. ([5]) in leading 1/iV-order, 
but it also contains corrections on the order of go/VN 
like Eq. pip . Such corrections, however, are irrelevant 
as they become only sizable at the Ehrenfest time te~ 
the time beyond which the classical solution fails and 
the true time-dynamics must be described anyway by the 
Schrodinger equation. 

B. Initial spin excitations 

Up to now we have focussed on a particular initial 
condition with excitations being initially present only in 
the boson mode. In contrast, a different initial condi- 
tion was considered in Ref. Il6l whereby the dynam- 
ics starts from an initial state with no boson present 
but, say, with the i th spin excited. Considering ho- 



mogeneous systems, it was found that during the time- 
evolution this i th spin remains excited if the total num- 
ber of spins is large, regardless of how strong the spin- 
boson coupling is. The corresponding expectation value 
is (0, lfti\S* (*)|0,4ti> = 1/2- (1 - cos (Qt)) /N. This re- 
sult was associated with the effect of "radiation trapping", 
and it can also be obtained with the classical approxi- 
mation. The corresponding initial condition Cj^i (0) = 
(0, 0,-1/2), d (0) = (0, 0, 1/2), and a (0) = is a fixed 
point of the classical Eqs. (0 [5]). Indeed, the effective 
magnetic field for each spin Bj = (0, 0, ej — ui) has only 
a z-component, and therefore the vector-product of two 
parallel vectors vanishes, i.e. B 3 x Cj = 0. The dynam- 
ics of the classical field a is frozen as Cj — 0. Quan- 
tum corrections to this result show up only in the first 
1/iV-correction. Thus, the classical approximation is also 
valid for a different initial condition in the few-excitation 
regime. 



VIII. CONCLUSIONS 

In this paper we have shown that the solution to the 
classical Hamilton equations of motion of the inhomo- 
geneous Dicke model coincide with the solution of the 
time-dependent Schrodinger equation when the number 
of spins N is large and the number of excitations p is 
small, p <C N. For a single excitation the leading 1/N- 
order of the classical solution coincides with the quan- 
tum solution. For a few excitations such correspondence 
does not hold but for p = 2,3 excitations the numeri- 
cal solutions of both classical equations of motion and 
Schrodinger equation coincide when the number of spins 
is large. It is plausible to conjecture that the same cor- 
respondence holds for p > 3 in leading 1/iV-order. 

To assess the validity of the classical approximation for 
p > 1 excitations, we have solved the Schrodinger equa- 
tion exactly in the two-excitation subspace with inhomo- 
geneous coupling constants only and compared the result 
with the classical solutions. For large N, we performed an 
1 /iV-expansion of the solution to the Schrodinger equa- 
tion and recovered the classical solution in leading order. 
Subleading 1/iV-corrections cause small deviations of the 
classical from the quantum solution, that, at a large time 
scale, make the difference between the two significant. 
This defines the Ehrenfest time that we identify in the 
limit of p < N as t e = y/N/ (g 2 ). 

Analyzing the solution to the Schrodinger equation 
for p — 2, we compared it with the solutions to the 
Schrodinger equation for p = 1. We have found that 
the boson occupation number in the two-excitation sub- 
space exhibits a multi-frequency dynamics due to the in- 
homogeneous couplings only, which, unlike in the single- 
excitation subspace, can lead to a decay in the limit of 
large N. But the leading term of an l/A^-expansion re- 
covers the single-frequency dynamics. The decay due to 
the inhomogeneity shows up only in the first sublead- 
ing 1 /TV-correction. We find that this contribution is an 
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oscillatory mode with frequency |17 and a slow decay en- 
velope. The decay is essentially non-exponential with a 
long power-law tail, and the decay time is r g ~ \/N/ (g 2 ), 
where W (g 2 ) is a characteristic coupling, the numerical 
prefactor is of order one for the special case of uniformly 
distributed coupling constants. 

The decay due to an inhomogeneous coupling to a spin 
bath, which is unavoidable as the spins are located at 
different positions of the cavity mode (with different am- 
plitudes of the electromagnetic field), is similar to the 
decay of an electron spin coupled to a bath of nuclear 
spins through the hyperfine interactio n 24 ' 25 . In the dy- 
namics of a cavity mode this mechanism can be neglected 
when only a few excitations are present in the system (for 
instance, in the few-photon spectroscopy experiments in 
Ref. [l|), but may lead to a significant decay in a sys- 
tem with many excitations present initially (such as, for 
instance, in a bath of nuclear spins coupled to a cavity). 



Expanding the above expression in 1/N, we obtain the 
corrections da (0) and 5(3 k (0) in Eqs. ([2"2" 1 12"5 ]| . 

To calculate the 1/iV-correction when £ = go, i.e. 
for maximum inhomogeneity, we expand Eq. (|19[) up 
to the second order in 1/N at the poles, s — ±iso, 
so = \/4N (<? 2 ), and also account for the second order 
corrections that come from the positions of the poles, 

so = ^N(g 2 )-lg 2 . 

Performing this procedure we write the residues of 
a (s) at s = ±iso as 

Res s=±4So a (.) e- = J_^±*£ e ±«M (A5 ) 



where 



A> = - 



±is D + SD 



2 (±igg) Ng 2 
(s 2 + N(g 2 )) 



(A6) 
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N = l + 



Ngl 



3(-s 2 + N{g 2 }) 



(A7) 



are the denominator and numerator obtained using just 
so = y/AN (g 2 ). Further, 



Appendix A: 1/N corrections to a p (t) and (t) 

In this appendix we calculate the first 1 /N correction 
to the pole contributions to the inverse Laplace transform 
of a (t) and f3 k (t) from Sec. V. 

If £ = 0, i.e. for homogeneous coupling constants, the 
solution to Eq. (I13[) simplifies. Substituting gj = go into 
Eqs. (dam, we get 



and 



M«) 



a(s) 



i\/2g 



+ (AN-2)g 2 -' 



ls 2 + (2N-2)g 2 



(Al) 



(A2) 



s s 2 + {AN -2)gl' 

These expressions have only poles, but no branch points: 
Two poles for /3k (s), s = ±i2go\/ N — |, and three for 



a(s), s — 0, H2go J N — |, see Fig. [2] The inverse 

Laplace transform is given by residues only, a (t) = a p (t) 
and (3 k (t) = 0> (t), 



isinfogoJN- It) 

ft (t) = V - V > (A3) 

k V ; V2A- 1 V ; 



and 



a p (t) 



+ — — — cos [ 2g \ N - -t | . (A4) 

2N- 1 2N- 1 \ y V 2 } K ' 



SD = - 



2(±is )Ng 2 



( _ s 2 +7V(5 2 )) 2 5(- S 2 +iV( 5 2 ))' 



SN 



2Ngl 



5(-s 2 + N(g 2 }y 



(A8) 



(A9) 



are the first 1/A-corrections. The average (g 2 ) = g 2 /3 
is evaluated using Qi (g) with £ = go- 
Performing the summation over the two poles Hso we 
get 



^ Res s=±iso a(s) 



cos(s t), (A10) 



±is 



?st = 2 N + 5N 
iso Do + SD 

and expand it in the small parameter as 
^ Res s=±iso a (s) e st 

±is 

2 N ( 6N\ ( SD\ , . /A . 
= —-^(l + — )(l- — )co S ( S ot), All 
isq D Q \ N Q ) \ Do) 

where the first two terms in the product still have to be 
expanded in the small correction to sg — W4N (g 2 ), and 
the last two have to be calculated using only the leading 
term s = \/4:N (g 2 ). 

The first two terms, which have to be calculated with 



,s = yj4N(g*)-§gl 



13 



_2_iVo 
iso D 



if J_ 

2\ + 107V 



20N (g 2 
6 



-9lN 



1 - 



5N 



2 (1 - 



2 

57V / i2y/I(g^ 
3 



1 - 



57V 



207V 



107V 



(A12) 



and the last two, which have to be calculated with sq = 
y/AN (g 2 ), are 

(' + £)('-£)-M)M)-4 

( Al3 ) 

Finally, the contribution from the poles ±iso is 

Res s=±J ; So ae st = - cos ( s ot)+^^ cos (s i) • (A14) 

±is 

The residue of a (s) at s = is also expanded in the 
small corrections, 



7V + ()7V 7V / <57V\ / 573 
Res ^ oa (s) = D^TW D~o\ TVo ) I A) 



where 



57V _ 9 5D _ 36 
TVo ~ 57V' Aj ~ 57V' 



(A15) 
(A16) 



And the contribution from the pole s — is 

Res s=0 a (s) 



2 207V 



(A17) 



The sum of Eqs. (|A14| IA17|) gives the correction 6a (go) 
from Eq. (|2"2j). 

Then, we calculate the corrections to /3jP (t) = 
Res,s = ±j S0 /3fc (s) e st . Expanding the residues at s — ±zs , 



(A18) 



where s = \/47V (g 2 ) - we get 



Res s=±iso /3 fc (s) e s 



-i\/2. 9fc e ±is °* 1 
s% ~ N (g 2 ) + 2g 2 D + 6D 



-iV2g k e ±lSot 1 ( 5D 



s 2 -N(g 2 )+2g 2 D Q ^ D J ' 



where Dq and 6D have already been calculated, see Eqs. 
(|A6l IA8[) . Similarly to the calculation of a (s) we again 
expand and get 



Res„ ± „ ft (.) «- = ^ ( 1 " ^£ + 4) , "i — + (A19) 



tv^ ^ a 5Ay (±i)4y7v(py v sa; v sa; 

The sum over these residues, 

ovt„ \ i9kV2 ( ^ 2(g k /g f i 6 ^ ^ | 9 



* w - i7^?y I 1 - l 1+ 5Fj am <«*> 

~ ( 1 2(ff fc /g ) 2 3\ 

1 + — sin (s t) , (A20) 



V2NW) V ^ 

r 



is the correction from Eq. ([23)1 for £ = go. 
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